%load('biannualdata.mat');
data=biannualdata(biannualdata.cycle=='1',:);
datafac=table(data.fac_id,data.period,data.suma,data.sumemi);
datafac.Properties.VariableNames = {'fac' 'period' 'alloc' 'emi'};
datafac=unique(datafac,'rows');


period=unique(datafac.period);
period=period(16:25);
nperiod=length(period);
datafac=datafac(ismember(datafac.period,period),:);

facid=unique(datafac.fac);
nfac=length(facid);

%take those fac with all last 10 observations
samplefac=[];
for i=1:nfac
    testfac=datafac.fac(datafac.fac==facid(i));
    len=length(testfac);
    if len==(25-16+1)
        samplefac=[samplefac;facid(i)];
    end
end

datafac=datafac(ismember(datafac.fac,samplefac),:);

rawdata=datafac;



[datafac,supply_yearly]=artfac(rawdata);
[datafac_norm1,norm1_average]=data_norm1(datafac,supply_yearly);

%added this
facid=unique(datafac.fac);
nfac=length(facid);

%grouping
round_norm1_average=ceil(norm1_average*coarseness)/coarseness; %creating groups by rounding k to the nearest 1/50
[groupaverage,ika, ikc]=unique(round_norm1_average);% groupaverage = round_norm1_average(ika) and round_norm1_average = groupaverage(ikc).
%[facaverage,ifa,ifc]=unique(norm1_average); % facaverage=norm1_average(ifa), from 570 to 57;norm1_average=facaverage(ifc), from 57 to 570
[fac,ia2, ic2]=unique(datafac.fac); % fac = datafac.fac(ia2) and datafac.fac = fac(ic2).

%obtaining d and k
%fac_updatek=norm1_average;
fac_d=norm1_average;
ngroup=length(groupaverage)
for it=1:ngroup
    allq=datafac_norm1.q(ismember(round_norm1_average,groupaverage(it)));
    fac_d(ismember(round_norm1_average,groupaverage(it)))=min(allq); 
    %fac_d(ismember(round_norm1_average,groupaverage(it)))=prctile(allq,minprct); 
end
datafac_norm1.facd=fac_d;
mini_demands=fac_d(ia2);
reserve=sum(mini_demands);
norm2=1/(1-reserve);
datafac.q=(datafac_norm1.q-fac_d)*norm2; %percentage of the allocable permit,sum to 1
datafac.d=fac_d*norm2; %minimum demand
mini_demands=datafac.d(ia2);
reserve=sum(mini_demands);
datafac.hatR=(1+reserve)*ones(size(fac_d));
datafac.r=norm2*datafac_norm1.r;
fac_updatek=norm1_average;
fac_ea=norm1_average;
for it=1:ngroup
    allq=datafac.q(ismember(round_norm1_average,groupaverage(it)));
    fac_updatek(ismember(round_norm1_average,groupaverage(it)))=prctile(allq,maxprct); %take the 90th percentile of all demand 
    fac_ea(ismember(round_norm1_average,groupaverage(it)))=mean(allq);
end
fac_indivk=zeros(nfac,1);
for indfac=1:nfac
    ref=fac(indfac);
    allq=datafac.q(datafac.fac==ref);
    maxq=prctile(allq,maxprct);
    fac_indivk(indfac)=maxq; %take the 90th percentile of all demand 
end
group=unique(fac_ea);
size(group)
groupk=fac_updatek(ika);
updatek=groupk;
grpsize=zeros(ngroup,1);
ea=zeros(ngroup,1);
for g=1:ngroup
    grp=group(g);
    grpsize(g)=length(datafac.fac(ismember(round_norm1_average,groupaverage(g))))/10;
    ea(g)=grp;
end


samplesize=length(unique(datafac.fac));

